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Abstract: We introduce a numerical procedure to investigate the spectrum of massive 
modes and its contribution for gravity localization on thick branes. After considering 
a model with an analytically known Schroedinger potential, we present the method and 
discuss its applicability. With this procedure we can study several models even when 
the Schroedinger potential is not known analytically. We discuss both the occurrence of 
localization of gravity and the correction to the Newtonian potential given by the massive 
modes. 
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1. Introduction 

In this work we deal with (4+l)-dimensional actions describing a scalar field coupled to 
gravity which leads to brane solutions with a sing le extra dimension 0, |, |, |, |, |]. As 
one knows, after an Ansatz is given for the metric, minimization of the action leads to a 
system of second-order differential equations for the field and the metric parameter related 
to the warp factor. With a convenient choice of the potential one can obtain a system of 
first-order equations, which helps to find analytic solutions for the metric and field. 

In this scenario, fluctuation around the solutions for the metric and scalar field can 
be decoupled in the transverse traceless gauge. After dropping the (3+l)-dimensional 
plane wave components, considered to satisfy a Klein-Gordon equation, the extra-dimension 
component of the metric fluctuations is reduced to a Schroedinger-like equation with an 
infinite number of solutions. Each solution is assigned to a massive mode characterized 
by the mass-shell condition imposed after the separation of the contribution of the four 
standard dimensions from the fifth or extra dimension. 

In order to obtain the Schroedinger-like equation one usually goes to a conformal 
variable related to the extra dimension. However, it occurs that such a transformation 
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is not always possible to be obtained analytically, and this may lead to a Schroedinger 
potential which is not known analytically for most of problems considered in the literature. 
Also, even when one knows the Schroedinger potential analytically, the whole spectrum of 
massive modes is hard to be obtained and normalized properly in order to study the main 
correction to the Newtonian potential. 

In the literature, Ref. J|| has already shown that any of the following two conditions 
is sufficient to prove the occurrence of gravity localization, namely one has to: i) consider 
the zero-mode, its normalization and asymptotic behaviour or ii) analyze the asymptotic 
behaviour of the Schroedinger potential. For instance, both procedures were used in Ref. 
Q to confirm gravity localization in a specific class of models. 

However, besides confirming the occurrence of gravity localization in a model, it is 
of physical significance to estimate the first-order corrections to the Newtonian potential. 
This can be related to the limit imposed by experiments on gravity at short distances ||. 
In this direction, Ref. j|] also presented an asymptotic analysis for the massive modes 
in a class of potentials that falls off as a(a + l)/z 2 far from the brane; as a result, a 
correction of the order of l/i? 2a_1 is found. The issue is that corrections to the Newtonian 
potential for analytical Schroedinger-like potentials can in principle be found in a similar 
way. However, it is not rare to face problems in which the Schroedinger-like potentials are 
only known numerically, and in such cases the above procedure cannot be used anymore. 
This poses the question on how one can generate and normalize properly the massive modes 
in order to obtain the corrections to the Newtonian potential. The issue has motivated us 
to investigate the problem, and to propose the present procedure. Our numerical study 
shows that in order for the spectrum of massive modes to be obtained properly one must 
apply the normalization procedure with great care. A wrong choice of the plane wave 
normalization can lead to a misunderstanding of the pattern of the distribution of the 
massive modes and to a wrong prediction of the corrections. 

In the present work, in the next Sec. || we deal with the equations one needs to build 
a Minkowski brane and to study the corresponding fluctuations in the metric and scalar 
fields. In Sec. | we review the analytic model of Ref. || and there we introduce the nu- 
merical method for the corrections of the Newtonian potential. A comparison of our result 
and the asymptotic analysis of Q is done in order to evaluate the precision of the method. 
We then discuss in Sec. [| some models with analytically known Schroedinger potential 
where we apply the known asymptotic analysis method, and in Sec. [5] we study models 
with numerically known Schroedinger potential, where we apply our numerical method. A 
comparison between the results coming from the numerical method and the asymptotic 
analysis is done. Also a comparison between the results coming from the numerical Nu- 
merov method and a simpler Euler method also available for non-analytic potentials is 
done. We end this work in Sec. [| where we include some conclusions and perspectives of 
future investigations. A short review of the Numerov method for determining the massive 
modes is presented in an Appendix. 
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2. The brane equations and fluctuations 



We start with 

S = J d^xdyy^^R+^d^-V^)] , (2.1) 

where g = det(g a b) and the metric 

ds 2 = g ab dx a dx b = e 2A r]^ u dx^dx v - dy 2 (2.2) 

describes a background with 4-dimensional Poincare symmetry with y as the extra dimen- 
sion. Here a, b = 0, 1,2,3,4, and e 2A is the warp factor. We suppose that the scalar field 
and the warp factor only depend on the extra coordinate y. 

The action given by Eq. Q2.1| ) leads to the following equations for the scalar field 4>{y) 
and the function A{y) from the warp factor: 

" + 4^' = ^M ( 2 .3a) 

dcp 

A" = -\cjj 2 (2.3b) 

A 12 = V 2 - \V{4>) (2.3c) 
o 3 

where the prime is used to represent derivative with respect to y. 
The potential is supposed to have the form 

where W(4>) is in principle an arbitrary function of the field eft - in the supersymmetric 
context W is named superpotential. The particular relation between V and W in ( $2.4 ) 
leads to a description in terms of a set of first-order differential equations, which are given 

by 1 1 

/ 1 dW . 
*=2l* < 2 ' 5a) 

A' = --W (2.5b) 

We consider now the effective 4-dimensional gravitational fluctuations in the confor- 
mally flat background discussed previously, as well as the fluctuation of scalar fields around 
solutions of the first-order equations (|2.5|) ; that is, we write: 

ds 2 = e 2A{y) {rj/u, + eh PbV )dx^dx u - dy 2 (2.6) 

and we set 4> — ► 4> + e< P where h^ u = h^ u (x,y) and </> = (f)(x,y) represent the fluctuations 
and e is a small parameter. On the transverse traceless gauge the metric perturbation 
separates from the scalars ||, leading to 

/^ + 4A'/^ = e - 2A DV, (2.7) 
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with □ being the (3 + l)-dimensional d'Alembertian. 

This equation can be decoupled by separating the 4-dimensional plane wave pertur- 
bations from the extra dimension contribution. We introduce a new variable z that turns 
the metric into a conformal one. This changes the equation for the extra dimension contri- 
bution of the metric perturbations in a Schroedinger-like form, where no single derivative 
terms are present. The new conformal coordinate z is defined by 

dz = e- A{y) dy. (2.8) 

The separation of variables is taken as 

h^(x,z)=e^ x e-l A ^^(z), (2.9) 

which turns Eq. ( |2.7| ) into a Klein-Gordon equation for the 4-dimensional components of 
the transverse-traceless Lj,, with the remaining Schroedinger-like equation 

" + Vsch(z)^ m (z) = m 2 4, m (z) (2.10) 

where we have dropped the [iv indices from the wavefunctions, now labeled by the corre- 
sponding energy m 2 . Here the Schroedinger-like potential is given by 

V sch (z) = ~A»(z) + ^A' 2 (z) (2.11) 
and this ends our revision of the standard braneworld scenario. 

3. An analytic volcano potential 

Before going into the issue concerning our numerical procedure, let us first consider the 
construction of a simple analytic potential characteristic of gravity localization. According 
to for potentials constructed with the procedure of Eq. ( |2.11| ) that go to infinity as 1/z 2 
have this property. Also, in order to have a volcano structure we need the two components 



of the rhs of ( 2.11 ) with opposite signs. A simple choice for A{z) is 

A(z) = -ln(l + z 2 ), (3.1) 
which is depicted in Fig. [|. In this case Eq. (grp 

gives the related Schroedinger-like 

potential 

V 3 , (OO) 

~~ + + (1 + ^)2" V-V 

which is also shown in Fig. [I] 

The characteristic volcano shape of the potential leads to the thick brane solution, 
so it contributes to smooth the potential considered in the original Randall-Sundrum Q 
model which describes the thin brane solution. The analytic structure of V sc h(z) eases the 
investigation of gravity localization. In particular, the zero-mode is determined analitically 
as 

N 

MZ) = (1+^)3/2 ' ( 3 ' 3 ) 

where Nq is a normalization constant. 
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Figure 1: Plot of A(z) (left) and the Schroedinger-like potential V sc h(z) (right). 
3.1 Reconstructing the model 

Before turning to the gravity localization issue, let us first reverse the investigation and 
work on the reconstruction of the model. This means to obtain the functions A(y), z(y), 
4>{y) and V((f>) once A(z) is given. We hope that this reconstruction process will give us 
insights on the importance of V sc h(z) in the far-reaching region for gravity localization. 
We first consider Eq. and write 

e~ A(y) dy. (3.4) 

Suppose we can write analytically z(y) in order to change the representation as A{z) = 
A(z(y)) = A{y). In this way we can invert Eq. fl2,8|) and write, in the present case 

z A ^dz = j e- h ^ 1+z ^dz = tan- 1 (z). (3.5) 

We invert this relation to write 

z = tan(y) (3-6) 



which we plot in Fig. We use it in Eq. (3.1) to get 

A(y) = -ln(l+tan 2 (y)). (3.7) 

which shows - see Fig. || - that the y variable is allowed to vary continuously only when 
spanning finite intervals, characterizing an effective compactification of the extra dimension 
in this representation. 

The warp factor is shown in Fig. ||[ There we see that it goes to zero far from the 
brane, as usual, but here this happens in a finite distance y* from the brane, achieving the 
characteristic AdS$ space in this region. 

The formation of a natural barrier on the extra dimension y must also be seen in the 
solution of the scalar field 4>{y). In order to recover this, we turn attention to Eqs. ([2.5]). 
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Figure 2: Plots of the functions y(z) (left) and A(y) (right). Note the effective compactification 
on the extra dimension y. 




Differentiating Eq. Q3.7| ) with respect to y and substituting in Eq. ( 2.5b| ) gives W{4>) = 
6tan(y). We can differentiate this to get W'(<p)4>'(y) = 6sec 2 y. We now use Eq. ( |2.5a| ) to 
obtain 

<f)(y) = \/31n(secy + tany). (3-8) 

In Fig. || one sees that there are finite intervals of y in which the function (j)(y) is not 
defined, and the brane has a diffuse structure. This can be better observed if we turn back 
to the z coordinate; we write 



<j){z) = y^) In (Vl + z 2 + z), (3.9) 

which is plotted in Fig. showing the kinklike profile corresponding to a difuse wall, as it 
appears in the vacuumless potential investigated in Q. 




Figure 4: Plots of <j>(z) (left), showing resemblance with a kink struture, but without an asymptotic 
value for 4>(z), and of z 2 V sc h{z) (right), showing the asymptotic regime for the Schroedingcr-like 
potential. 

In this case the superpotential can also be analytically determined in the following way. 
An ordinary differential equation of first order for W{4>) is obtained after using Eqs. (|2.5| ) 
and derivatives of A(y) in Eq. Q3.7|) and of <f>{y) in Eq. ( |3.8| ). The procedure leads to 
W'(cp) + l/V3W(4>) = 2\/3expc/>/v / 3. We solve this equation to obtain 

W(4>) = 3e 2/3< ^ - A (3.10) 

with A being a free parameter. The potential is then given by 

V{<t>) = --e^ + 2Ael^--A 2 (3.11) 
2 3 

The analytic model is very nice, and it leads us to go forward or backward very nicely. 
3.2 Massive modes 

Now we turn back to the Schrodinger-like equation with the aim to study the massive 
modes. As one knows, in the present scenario it is important to find the full spectrum of 
massive modes. The issue is of greater importance within the context of local localization 
of gravity [fR]]. Although it is not always possible to find analytical expressions for the 



massive modes |11], the analytical form of the potential V sc h helps us very much to control 
the numerical investigation. Indeed, a simple Runge-Kutta routine can be used with the 
boundary conditions ij} m (0) = 1 and ip' m (0) = 0. The first condition is arbitrarily fixed, and 
will be adjusted with the normalization procedure. As one knows, the conformal z— variable 
leads to a probabilistic interpretation for |^ m (z)| 2 and one must impose the condition 

\i> m (z)\ 2 dz = 1. (3.12) 

We normalize the wavefunctions in a box of size [— z max , z max ]. However, an important 
point is to consider a box with sufficiently large size z max in order for the potential V sc h 
to achieve the 1/z 2 regime. In fact, as already noted in Q, localization of gravity is 
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Figure 5: Plots of |V'm(0)| 2 as a function of m. Here each ip{z) has been solved by Maple with 
Runge-Kutta default method of fourth order, taking z max = 30 (upper), z max — 10 (lower, left) 
and z max = 2 (lower, right). All plots show the asymptotic limit corresponding to l/z max . 

determined by the far region of the potential. As we can see in Fig. |||, for the choice of 
V sc h the asymptotic regime is achieved for z > 20. In this way we must choose an at least 
comparable size z max for our box in order to correctly reproduce the gravitational features. 
However, to illustrate possible mistakes within the procedure, we will consider boxes of 
different sizes, to see spurious effects which can be observed for small boxes. We plot some 
results in Fig. [5], where we show the normalized |?/> m (0)| 2 as a function of m for z max = 30, 
Zmax = 10 and z max = 2. Note that for z max = 30 we achieved the asymptotic regime 
for the 1/z 2 behaviour for the V sc fo potential. On the other hand, for smaller values of 
Zmax, say z max = 10, besides a first small peak, small perturbations resembling resonances 
start to appear. This behaviour is greatly reinforced for z m ax = 2. Clearly there are no 
resonances at all in this problem, the effect just observed being a consequence of a wrong 
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choice of the size of the box, which could not accommodate the correct behavior of the 
potential V sc h- 

Let us now turn to asymptotic analysis in order to evaluate the degree of accuracy 
of our Fig. H|, for z max = 30. First of all, note that in this region the factor |^ m (0)| 
increases monotonically with m, achieving a plateu for larger masses. The very existence 
of this plateau is confirmed after considering to 2 >> V sc h in the Schroedinger equation, 



Eq. ( 2.10 ). Indeed, in this case it reduces to 

d 2 ip m (z) 



dz 2 



m ipm(z), m » V sch 



(3.13) 



with solution ip m (z) = Ncos(mz). After normalization in a box one obtains IV'mlO)] 2 = 
l/Zmax for large m, which for z max = 30 gives ip m (z) ~ 1/30. See Fig. |5| for the accuracy 
of this limit. 

Now let us consider the regime of lower modes. According to Ref. Q - see particularly 
Sec. 4 -, this regime is related to the z » 1 behaviour of V sc h- Indeed, in our case we 
have the limit V sc h ~ 12/ z 2 . This has the particular form a(a + l)/z 2 , already proposed 
in Ref. There, it is shown that when this occurs one expects the lower massive modes 
to obey the relation tp m (0) ~ m a ~ l . For our potential we have a = 3, and this leads to 
IV'mlO)! 2 ~ m 4 . In order to look for this power law behavior in |?/> m (0)| 2 , we plot in Fig. || 
the former results of Fig. |B| for z max = 30 in a logarithmic scale, focusing on the lower 
modes. 




Figure 6: Plot of log(|-i/'»n(0)| 2 ) as a function of log(m) (left panel). Each ip(z) was generated with 
500 points, taking z max = 30. We are also displaying the straight line corresponding to the limit 
|VVt(0)| 2 ~ m A . The right panel shows the same plot, after dropping out the data which come from 
very small values of m, since they come from the plane wave normalization in a box with size not 
large enough. We are also displaying the straight lines corresponding to the lower (green line) and 
higher (yellow line) masses. 

Note that the desired behaviour for small masses is obeyed in the model for almost two 
orders of magnitude. However, for very small masses, say for m < to*, the approximation 
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breaks down, indicating the formation of a first small peak already noted in Fig. || Our 
simulations show that the region of validity of the approximation |^ m (0)| 2 ~ tb 4 is enlarged 
toward smaller masses when the parameter z max is increased. Qualitatively, this can be 
understood very easily, since the wavefunctions have larger wavelenghts for smaller masses. 
As a consequence, the wavefunctions in this region do not oscillate sufficiently in the box 
and the procedure of normalization cannot be justified. With this reasoning, the procedure 
of carrying out the integration of massive modes to obtain the Newtonian potential must be 
done with great care, in order to eliminate the wrong contribution of modes for < m < m* . 
One possible solution, justifiable by our asymptotic analysis, is to extrapolate the observed 
polynomial behaviour toward the region < m < m*, as shown in the right panel of Fig. 0. 

3.3 Gravity localization: numerical procedure and asymptotic analysis 

Let us now investigate the Newtonian potential for two masses mi and 777-2 separated by a 
distance R. In order to obtain the potential one sums the tower of Kaluza-Klein excitations 
to the usual contribution coming from the zero-mode as § 

U[R) = G ^ + ^J~ dm £^ Hmm >, (3.14) 

where G = M^ 2 represents the four dimensional coupling, M* is the fundamental five 
dimensional Planck scale and the integration is done considering the brane position at 
z = 0, as it is done in the thin brane case. Separating the contributions of the Einstein- 
Hilbert action due to the four dimensional part and the extra dimension, we are led to an 
expression relating the two scales as Q 

Mf = M*Nl (3.15) 

where Nq = after the normalization procedure given by Eq. (|3.12| ) is applied to the 
zero-mode Eq. fl3.3| ). In this way, Eq. (3.14) can be written as 



3-7T r°° 

1 + — / dme- mR \^ m {<d)\< 
° Jo 



(3.16) 



Note that in order to obtain the Newtonian potential U for a given separation R, it 
is sufficient for practical purposes that the integration goes over the massive modes until 
masses of O(10/R). For higher distances R, a considerable simplification is achieved, but 
for lower distances one must consider higher massive modes. In the limit, for R — > we 
must integrate all possible masses until m — > 00. In this way, another extrapolation of the 
obtained data for |^ m (0)| 2 is necessary for higher modes. Namely, one needs to extrapolate 
the results displayed in Fig. ||, right panel, as |V> m (0)| 2 ~ l/z ma x- 

Let us now deal with gravity localization from the asymptotic analysis of | , m (O)| 2 . We 
can separate the regime |V'm(0)| 2 = C\m 2 ^ a ~ l \ valid for smaller masses from the regime 
|^ m (0)| 2 = C 2 /z 

max-, valid for larger masses. The constants C\ and C2 can be learned 
from Fig. |(| right panel. This figure shows that the transition region can be located at 
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log(m/r) ~ 1. In this way we can separate both contributions for the Newtonian potential 
as follows 



' " •"' ' dmt — 



1 + r dme-™ R m^ + *L f 

° J0 ° Jrr 



For the model under investigation we have a = 3, and this leads to 



rriT Zmax 



(3.17) 



U(R) * 



1 _ ^ ( -24 + e —*r4i? 4 + 4m|i? 3 e- w ^ 
8 \ 



+ Ylm\R 2 eT mTR + 24m T i?e" mTi? + 24e~ mT ^ + 



3^ 1 e~ mTR 
R 



and for mxR >> 1 we can write 



i + 3 V 24 



, (3.18) 



(3.19) 



which reproduces the Newtonian potential for large values of R. We then learn from this ex- 
pression that the present model reproduces the known gravitational limit at large distances, 
localizing gravity with a 1/i? 6 correction. 

However, the full integration of Eq. (|3.17| ) is still necessary in order to estimate the 
accuracy of our procedure. This is important in case the potential V sc h is not analytically 
know and also to find the value of the parameter C\. To implement this, we fit numerically 
the logarithmic plot of Fig. ||, right panel; we did that considering that it comes approx- 
imately as a superposition of two straight lines, one for each region where m < rriT and 
m > tut- We remark that this is important in order to better achieve a control for models 
where the Schroedinger potential are not known analytically. First of all we get from a 
linear fitting of the Fig. ||, right panel, the result log|V'm(0)| 2 = —7.186 + 4.359 log m, 
leading to |^ m (0)| 2 = 7.57 x lCT 4 m 4 - 359 , giving C\ = 7.57 x 10 -4 . From the same figure 
we can estimate mx as the value of m corresponding to the crossing of both straight lines 
for larger and smaller masses. This gives log(m^) = 0.868, or my ~ 2.383. Also, we find 
for the region m > hit the limit log | t/j m (0)\ 2 = l/z max ~ 1/30. After substituting these 
results in Eq. ( |3.17| ) and performing the integration, we can find the correction coefficient 
L for U — U Newton ~ 1/R L , which is displayed in Fig. [7|. The numerical results give a 
correction L ~ 6.34 for large distances, which is to be compared with L = 6, as the value 
which comes from the asymptotic analysis. 

3.4 The Numerov method 

We can solve the same problem for the massive modes using the Numerov method |13[| 
instead of the Runge-Kutta fourth-order method that we have just used. The Numerov 
method - see Appendix - is largely used in order to look for bound states in quantum 
mechanics. We can also apply it here, and some results for this model are displayed in 
Fig. |8[ Despise the slow convergence of the Numerov method compared with Runge-Kutta 
methods, we can use the Numerov method when the Schroedinger potential is only known 
numerically. 
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The right panel in Fig. ^ shows a good 
visual concordance between both methods 
for |-0 m (O)| 2 as a function of m. However, 
in the former subsection we have empha- 
sized the importance of the lower massive 
modes for gravity localization. A careful 
analysis in a logarithmic scale has shown 
that for small masses, the discrepancy be- 
tween the methods grows, with Numerov 



log(R) 



method tending to present |^ m (0)|" 



~ m 



P 



with a higher power law j3 than the Runge- 
Kutta method. Naturally, this may reflect 
on predicting a correction to the Newtonian 
potential with a higher power law. 

We will use this method in Sec. | for 
some other models, where the Schroedinger- 
like potential is only known numerically, since there the Runge-Kutta method cannot be 
used anymore. 



Figure 7: Plot of the L coefficient as a function 
of the distance R between two unit test masses. 
The results are obtained after performing an ex- 
pansion valid for R>> 1/tot with hit — 2.38. 





Figure 8: Plots of tp(z) for m — 1 and z max — 30 (left). The lines are displayed to achieve a better 
comparison between the default method from Maple (500 points-black line), for the analytically 
known V sc h an d the Numerov method (2f2-red, 490-green and 980 points-blue lines) applied for 
z > 0, better indicated for numerically known potentials V sc h- Note the slow convergence of the 
method. The right panel shows |V>m(0)| 2 as a function of m, comparing Runge-Kutta method (black) 
with Numerov method for ip(z) with 212 points (6z = 0.14, red line) and 980 points (6z = 0.0316, 
blue line). 



4. Models with analytic V sc h- asymptotic analysis 



If we compare with the presented numerical procedure, the asymptotic analysis Q has very 
good precision to estimate the power of the first correction to the Newtonian potential. 
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However, the study can only be done when one knows an analytical relation for y = y(z). 
To make this point clearer, let us illustrate this possibility with the following examples, 
recently considered in Ref. [O]. 

4.1 The case Wx{4>) = 3asin(60) 
For this model we have 

A(y) = -2 A log(g COS h(^ab 2 y)). (4.1) 

In general this model does not give an analytic relation for z(y). However, for b 2 = 2/3 we 
get 

z = f e - A{y) dy = -sinh(ay), (4.2) 
J a 

with 



' z 2 a 2 



A(z) = -log( g y-^--l). (4.3) 
This model gives an analytic expression for V sc h, depending on two parameters q and a: 



2 (5z 2 a 2 + 2q 2 ^ 
4(g 2 — a 2 z 2 



V sch = 3a 2 (4.4) 



In this way, we can apply the asymptotic analysis in order to study the Newtonian potential. 
In fact, we have the limit 

15 , . 

Vsch = 4^2 > z » L ( 4 - 5 ) 

We see from the former expression the independence of the potential of the parameters q 
and a for large values of z. Following the former section we can write V sc fo = a(a + l)/z 2 
for z >> 1. Thus, a = 3/2 and we expect that |V'm(0)| 2 ~ m for small masses. This gives 
the same correction 1/E? for the Newtonian potential as in the Randall-Sundrum model. 

4.2 The case W 2 (4>) = 3asinh(&</>) 
For this model we have 

^) = -^log(sec 2 (^)). (4.6) 

For b 2 = 1/3 we can obtain an analytic expression for z(y), namely 

z = -tan(^). (4.7) 
a 2 

Inverting this expression and substituting in the expression for A(y) we get 

A(z) = -log[l + ^f). (4.8) 
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For this model we have an analytic expression for the Schroedinger potential 



1 2a 2 (z 2 a 2 - l 
(z 2 a 2 + 4) 

with the limit 



Vsch = 7-2 2 . ^ ( 49 ) 



Vsch = ^>, z» 1. (4.10) 

z z 

This limit leads to the same asymptotic result of the former section, where a = 3 and 
I (0)| 2 ~ m 4 , with a 1/ii 6 correction to the Newtonian potential. 



5. The case of models with numerical V sc h 

For a large class of models there is no analytic solution for y = y(z). In this way, the analysis 
of A(z) and the potential V sc h must be done numerically with no a priori asymptotic 
analysis available to guide us. In this case the method presented in this work may be used 
to estimate the behavior of |^ m (0)| 2 for small masses and to find the first correction to the 
Newtonian potential. We illustrate the procedure with some examples. 

5.1 The case W 3 {f) = 2atan" 1 (sinh6</>) 
This model gives (l!| 

A(y) = -Llog(l + a 2 6V) - \ay twT^atfy). (5.1) 

We choose a = 1 and write c = b 2 to get 

A(y) = ^-log(l + c 2 y 2 ) - -ytajT\cy). (5.2) 
3c 3 

and now we study the influence of c on gravity localization. In this case, even the relation 
given by Eq. ( p.4[ ) for z(y) is not possible to be obtained analytically. Thus, we have in- 
tegrated it numerically for some values of c, leading to a numerical determination of A(z) 
which we show in Fig. ||. Gravity localization is confirmed after demonstrating normaliza- 
tion of the zero-mode, as we also show in Fig. |9[ We have repeated the same analysis in a 
logarithmic scale, and the results shown that the zero mode is nicely normalized as it goes 
to zero faster than z~ 1 / 2 , thus confirming the occurrence of gravity localization. 

In order to study details of gravity localization one must construct the potential V sc h- 
We did this in Fig. |l0|, extending the z axis in order to guarantee that the asymptotic regime 
is correctly achieved and that the higher massive modes can be properly normalized. Fig. |l0| 
shows that even without an analytic expression for V sc h, it indicates that asymptotically 
one has V sc h = a ( a + l)/ z2 > an d that this regime is better achieved for larger values of c. 
In particular, for the case c = 0.5 we have found for the limit of V sc hZ 2 the value 4.02(2) 
which gives a = 1.52(4). 

In Fig. 0, upper panel, we show the normalized probability on the brane 1^/^(0)1 as 
a function of m. We see that for lower values of c (red thin line) the influence of massive 
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Figure 9: Plots of A(z) (left) for the model with W±{(j)) = 2tan~ 1 (sinhfo</>), for the values c = 0.5 
(red, thin), 1 (green, thick) and 2 (blue, thicker), with Sz = 0.05. We also show the zero-mode 
tpo{z) ~ e 3A / 2 (right) as a function of z for the same values of c, compared with the z~ x / 2 curve 
(black). 




Figure 10: Plots of V sc h{ z ) (left) for W±{(f>) = 2a tan 1 (sinh&</>) and z 2 V sc h(z) (right), showing that 
the asymptotic regime is nicely achieved for larger values of z. We are using the same conventions 
of Fig. 9. 

modes increase faster than for higher values of c (blue thick line). This effect shows that 
gravity is easier localized for higher values of c, in agreement with the results shown in Fig. 

1 

Now, the behavior of V sc h(z) for z >> 1 leads to |VVra(0)| 2 ~ m 2 ( a_1 ) ~ m 1 ' 1 ^ 1 ' for lower 
masses. We compare this result with the one shown in Fig. [ll], left, the green line which 
represents the best fit for lower masses, obtained after neglecting the very low points. In 
this figure, we show the logarithmic plot of |^ m (0)| 2 as a function of m, after considering 
z max = 200 with step 5z = 0.05, using the Euler method for determining the massive 
modes. Note that, as the number of iteractions increased, it is important to find with 
higher precision the y points corresponding to z with fixed step. The plot indicates the 
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Figure 11: Plot of \ip m (0)\ 2 (m) (upper panel), using y 1 = V0.02 ~ 0.14 with 700 points (Sz = 0.14) 
leading to y m ax — 4.72, using the same conventions of Fig. 9. In the lower panel we plot refined 
analysis for c = 0.5, showing logarithmic plots of \ip m (0)\ 2 (m) for z max ~ 200 with 4000 points, 
using Euler method (left) and Numerov method (right), fixing the step on z with precision of 10 -7 . 



appearance of a power law |V>m(0)| 2 ~ mP for lower masses, with ^Euier = 1-3, or after 
using P = 2(a - 1), a Eu ier = 1-65. 

The very same simulation with Numerov method gives j3 'Numerov = 1.6 or a Numerov = 



1.80 - see Fig. 11, right panel. We notice that the Numerov method produces consistent 
results both for higher and lower masses, as we would expect from a higher precision 
method. On the contrary, the Euler method deviates from expected results in those limits 
due to error propagation. However, as the simpler Euler method gives fast results, it was 
used as a guide to search for the adequate regions to apply the Numerov method. In this 
way, we found two complementary numerical estimations for the correction 1/R L for the 
Newtonian potential for large distances, L = j3 + 2, leading us to the results L Euler ~ 3.3 
and LNumerov ~ 3.8. These results are to be compared to Lschroedinger ~ 3.1(1), which was 
produced by asymptotic analysis of the Schroedinger potential. 
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5.2 The case W 4 (0) = (l/2)m0 2 + b 
For this model we have 

My) 



l e my _ by 

6 3 ' 



In general one cannot obtain an analytic expression for z(y). However, for b 
use the exponential integral to write 

z = --Ei{l- e —) 
m o 



(5.3) 



we can 



(5.4) 




but we have been unable to invert this expression in order to get an analytic expression 
for y(z). 

In general, this class of models leads to 
an asymmetric warp factor. In particular, 
for negative values of b the warp factor goes 
to zero at a negative z value whose mod- 
ule increases for smaller values of b. An 
illustration of this is given in Fig. 12 for 
b = —1, in which one depicts the obtained 
warp factor. The figure shows the asymme- 
try of the constructed brane, and also the 
influence the m factor induces on the posi- 
tive z region, with a maximum warp factor 
that increases for lower values of m. The 
case of asymmetric brane has been studied 
in 0. 

The asymmetry of this model poses in- 
teresting questions concerning the influence of the asymmetry on gravity localization, since 
it leads to an asymmetric function z(y) as well as an asymmetric V sc h- Usually, the massive 
modes tend to acquire the characteristics of plane waves far from the brane and the nor- 
malization for symmetric models is realized after considering a symmetric box with length 
L. For asymmetric models, however, the normalization procedure for the massive modes 
must be done with care and will be postponed to a future investigation. 



Figure 12: Plots of e 2A(z *> for the W 3 (<f>) model, 
for the values b = — 1 and m — 0.5 (red, thin), 
m = 1 (green, thick) and m = 2 (blue, thicker). 



6. Conclusions 



In this work we have investigated gravity localization in diverse braneworld models. In 
general, we can classify the models into two distinct classes, one which leads to analytic 
Schroedinger-like potential, and the other, in which the Schroedinger-like potential cannot 
be given analytically. The first class of models is very important, and it was nicely studied in 
Ref. [|| , where an interesting procedure to fully investigate stability, with a very nice recipe 
to qualify gravity localization and Newtonian correction and to quantify the Newtonian 
correction. 
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However, there are interesting models which lead to Schroedinger-like potentials that 
cannot be given analytically. For such models, however, the method developed in 
cannot be used to obtain the massive models, thus inducing a gap on the study of gravity 
localization in braneworld models. This problem is yet more important within the context 
of local localization of gravity [10], since there the knowledge of the massive modes is crucial 



to understand the local localization of gravity. This fact has inspired us to propose the 
present study, to close the aforementioned gap in the braneworld proposal. 

In the course of our investigations, we had to face several interesting issues, among 
them we would like to pinpoint the following: 

1. In the Schroedinger-like equation, in order to correctly normalize the eigenfunctions 
we have to properly choose the size of the box, the z max . As we have shown, we have 
to deal with this with care, in order to prevent the inclusion of fake resonant states. 

2. When the Schroedinger-like potential may be obtained analytically, a comparison 
between the Runge-Kutta and Numerov methods was done. In this case, the Runge- 
Kutta method was applyed directly as a default method for Maple, with an absolute 
error that can be easily fixed to a small value. As expected, the best fit between both 
methods is for Numerov method with smaller steps 5z. 

3. When the Shroedinger-like potential is analytically known, the asymptotic analysis 
is the best choice in order to obtain the power law of the first correction to the 
Newtonian potential. 

4. Despide the Numerov method being the best method used to find bound states, its 
precision is reduced when one changes from a boundary condition problem to a initial 
value problem. However, its stability is best seen when compared to Euler method. 
Also, the Runge-Kutta method appears to be inadequate when the potential is not 
known analytically. 

We believe that the present study contributes to enlarge the scope of stability in the 
braneworld scenario in the presence of large extra dimensions. The importance of the 
numerical procedure widens with the current interest on asymmetric branes, since there 
the extra dimension drives the system to behave differently, depending on the positive or 
negative sense one spans the fifth dimension. We hope to return to this subject in another 
work, investigating the presence of massive modes for asymmetric branes. 
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Appendix: Construction of eigenfunctions of Schroedinger equation by 
Numerov method 

Consider the Schroedinger equation 

= [V sch (z) - m 2 ]i( Z ) (Al) 

Given an energy m 2 , we look for the corresponding function ip m (x), when the potential 
V sc h(z) is an even function known only numerically. We follow the theoretical analysis 
of the renormalized Numerov method |l3|], where the linear character and the absence 
of the first derivative on the Schroedinger equation conspire to achieve a powerful sixth- 
order numerical method. Then, defining f(z) = V sc h(z) — m 2 , we can rewrite Eq. (A.l) 
as tp"(z) = f(z)ip(z). We consider the potential V sc h(z) given as a set of points, with 
the variable z described with fixed step h. This is very important and in the present case 
where the relation y(z) - see Eq. ( |2.8| ) - cannot be described analytically one needs a 
specific numerical procedure, where a convenient variable step on y leads to a constant 
step on z. We then expand the functions ip(z + h) and ip(z — h) in Taylor series in order 
to get 

4>(z + h)+^{z-h) = 2^{z) + /i V 2) (z) + ^V (4) (z) + ^/iV (6) (*) + h.o.t. (A2) 

12 Sou 

On the other hand, the second derivative of both sides of the former equation (multiplied 
by the convenient numerical factor — (l/12)/i 2 ) leads to 

- ^h 2 W'{z + h)+ r(z - h)} = -\h 2 ^ 2 \z) - 1/*V (4) - ^/»V (6) (*) + h.o.t. (A3) 

Adding both Eqs.(A.2) and (A. 3) after substituting (A.l) on the terms where the second 
derivative appears, we get 

ip(z + h) + %j){z - h) - 

^h 2 [f(z + h)i>(z + h)+ f(z - h)^(z - h)] = 2^{z) + h 2 f(z)^(z) + h.o.t. (A4) 

We now define T(z) = (l/12)/i 2 /0) to obtain the recurrence formula for ip(z): 

[1 - T(z + h)]ip(z + h) + [l- T(z - h)]^{z -h) = 2[1 + 5T(z)}<<p(z) (A5) 

We consider that the boundary condition is ^(O) = 1 and that we are searching only 
for even modes ip m {z). Then, we apply the recurrence formula for z = 0: 

[1 - T{h)]i>(h) + [1 - T(-h)]if>(-h) = 2[1 + 5T(0)]V(0) (A6) 
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The condition that V sc h(z) is even gives T{—h) = T(h). In this way the former equation 
becomes 

1 + 5T(0) 

This result and the recurrence formula (A. 5) lead us to the following strategy: we divide 
the interval for z > in N slices, getting the points z±, Z2, ...zn- Defining ip(zi) = 1, we 
obtain if) fa) = (1 + 5Tfa))/(l — Tfa)). Thus, for 3 < i < N we construct the points 

^ Zi) = 2[1 + 5^-1)1^(^-1) - [1 - T(z,_ 2 )]^(z. t _ 2 ) 

1 — T(zi) 

in order to obtain the desired massive mode as a set of N points, as we use in the text. 
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